﻿/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * Geophysical Computational Tools & Library (GCTL)
 *
 * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
 *
 * GCTL is distributed under a dual licensing scheme. You can redistribute 
 * it and/or modify it under the terms of the GNU Lesser General Public 
 * License as published by the Free Software Foundation, either version 2 
 * of the License, or (at your option) any later version. You should have 
 * received a copy of the GNU Lesser General Public License along with this 
 * program. If not, see <http://www.gnu.org/licenses/>.
 * 
 * If the terms and conditions of the LGPL v.2. would prevent you from using 
 * the GCTL, please consider the option to obtain a commercial license for a 
 * fee. These licenses are offered by the GCTL's original author. As a rule, 
 * licenses are provided "as-is", unlimited in time for a one time fee. Please 
 * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
 * to include some description of your company and the realm of its activities. 
 * Also add information on how to contact you by electronic and paper mail.
 ******************************************************/

#include "regular_grid.h"

#ifdef GCTL_WAVELIB
// include wavelib headfile
#include "wavelib.h"
#endif // GCTL_WAVELIB

void gctl::regular_grid::init(std::string in_name, std::string in_info, int xnum, int ynum, 
	double xmin, double ymin, double dx, double dy)
{
	if (initialized)
	{
		std::string err_str = "The object is already initialized. From regular_grid::init(...)";
		throw runtime_error(err_str);
	}

	base_mesh::init(REGULAR_GRID, MESH_2D, in_name, in_info);

	rg_xnum = xnum; rg_ynum = ynum;
	rg_xmin = xmin; rg_ymin = ymin;
	rg_dx = dx; rg_dy = dy;

	node_num = rg_xnum * rg_ynum;
	ele_num  = (rg_xnum-1) * (rg_ynum-1);

	nodes.resize(node_num);
	elements.resize(ele_num);

	for (int j = 0; j < rg_ynum; j++)
	{
		for (int i = 0; i < rg_xnum; i++)
		{
			nodes[i+j*rg_xnum].id = i + j*rg_xnum;
			nodes[i+j*rg_xnum].x  = rg_xmin + i*rg_dx;
			nodes[i+j*rg_xnum].y  = rg_ymin + j*rg_dy;
		}
	}

	for (int j = 0; j < rg_ynum-1; j++)
	{
		for (int i = 0; i < rg_xnum-1; i++)
		{
			elements[i+j*(rg_xnum-1)].id = i+j*(rg_xnum-1);
			elements[i+j*(rg_xnum-1)].vert[0] = nodes.get(i+j*rg_xnum);
			elements[i+j*(rg_xnum-1)].vert[1] = nodes.get(i+1+j*rg_xnum);
			elements[i+j*(rg_xnum-1)].vert[2] = nodes.get(i+1+(j+1)*rg_xnum);
			elements[i+j*(rg_xnum-1)].vert[3] = nodes.get(i+(j+1)*rg_xnum);
			elements[i+j*(rg_xnum-1)].dl = nodes.get(i+j*rg_xnum);
			elements[i+j*(rg_xnum-1)].ur = nodes.get(i+1+(j+1)*rg_xnum);
		}
	}

	initialized = true;
	return;
}

void gctl::regular_grid::load_binary(std::string filename)
{
	if (initialized)
	{
		std::string err_str = "The object is already initialized. From void regular_grid::load_binary(...)";
		throw runtime_error(err_str);
	}

	std::ifstream infile;
	gctl::open_infile(infile, filename, ".2m", std::ios::in|std::ios::binary);

	// 读入网格头信息
	load_headinfo(infile, REGULAR_GRID, MESH_2D);

	// 读入网格信息
	infile.read((char*)&rg_xnum, sizeof(int));
	infile.read((char*)&rg_ynum, sizeof(int));
	infile.read((char*)&rg_xmin, sizeof(double));
	infile.read((char*)&rg_ymin, sizeof(double));
	infile.read((char*)&rg_dx, sizeof(double));
	infile.read((char*)&rg_dy, sizeof(double));

	node_num = rg_xnum * rg_ynum;
	ele_num  = (rg_xnum-1) * (rg_ynum-1);

	nodes.resize(node_num);
	elements.resize(ele_num);

	for (int j = 0; j < rg_ynum; j++)
	{
		for (int i = 0; i < rg_xnum; i++)
		{
			nodes[i+j*rg_xnum].id = i + j*rg_xnum;
			nodes[i+j*rg_xnum].x  = rg_xmin + i*rg_dx;
			nodes[i+j*rg_xnum].y  = rg_ymin + j*rg_dy;
		}
	}

	for (int j = 0; j < rg_ynum-1; j++)
	{
		for (int i = 0; i < rg_xnum-1; i++)
		{
			elements[i+j*(rg_xnum-1)].id = i+j*(rg_xnum-1);
			elements[i+j*(rg_xnum-1)].vert[0] = nodes.get(i+j*rg_xnum);
			elements[i+j*(rg_xnum-1)].vert[1] = nodes.get(i+1+j*rg_xnum);
			elements[i+j*(rg_xnum-1)].vert[2] = nodes.get(i+1+(j+1)*rg_xnum);
			elements[i+j*(rg_xnum-1)].vert[3] = nodes.get(i+(j+1)*rg_xnum);
			elements[i+j*(rg_xnum-1)].dl = nodes.get(i+j*rg_xnum);
			elements[i+j*(rg_xnum-1)].ur = nodes.get(i+1+(j+1)*rg_xnum);
		}
	}

	initialized = true;

	// 读入模型数据单元
	load_datablock(infile);

	infile.close();
	return;
}

void gctl::regular_grid::save_binary(std::string filename)
{
	if (!initialized)
	{
		std::string err_str = "The object is not initialized. From void regular_grid::save_binary(...)";
		throw runtime_error(err_str);
	}

	std::ofstream outfile;
	gctl::open_outfile(outfile, filename, ".2m", std::ios::out|std::ios::binary);

	// 首先输出网格的头信息
	save_headinfo(outfile);

	// 输出网格信息
	outfile.write((char*)&rg_xnum, sizeof(int));
	outfile.write((char*)&rg_ynum, sizeof(int));
	outfile.write((char*)&rg_xmin, sizeof(double));
	outfile.write((char*)&rg_ymin, sizeof(double));
	outfile.write((char*)&rg_dx, sizeof(double));
	outfile.write((char*)&rg_dy, sizeof(double));

	// 输出的模型数据单元
	save_datablock(outfile);

	outfile.close();
	return;
}

void *gctl::regular_grid::get_node(unsigned int idx)
{
	std::string err_str;
	if (!initialized)
	{
		err_str = "The object is not initialized. From gctl::regular_grid::get_node(...)";
		throw runtime_error(err_str);
	}

	if (idx >= node_num)
	{
		err_str = "Invalid index. From gctl::regular_grid::get_node(...)";
		throw runtime_error(err_str);
	}

	return nodes.get(idx);
}

void *gctl::regular_grid::get_element(unsigned int idx)
{
	std::string err_str;
	if (!initialized)
	{
		err_str = "The object is not initialized. From gctl::regular_grid::get_element(...)";
		throw runtime_error(err_str);
	}

	if (idx >= ele_num)
	{
		err_str = "Invalid index. From gctl::regular_grid::get_element(...)";
		throw runtime_error(err_str);
	}

	return elements.get(idx);
}

gctl::regular_grid::regular_grid() : base_mesh::base_mesh(){}

gctl::regular_grid::regular_grid(std::string in_name, std::string in_info, int xnum, int ynum, 
	double xmin, double ymin, double dx, double dy)
{
	init(in_name, in_info, xnum, ynum, xmin, ymin, dx, dy);
}

gctl::regular_grid::~regular_grid(){}

int gctl::regular_grid::get_xdim() const
{
	if (!initialized)
	{
		std::string err_str = "The object is not initialized. From void regular_grid::get_xdim(...)";
		throw runtime_error(err_str);
	}

	return rg_xnum;
}

int gctl::regular_grid::get_ydim() const
{
	if (!initialized)
	{
		std::string err_str = "The object is not initialized. From void regular_grid::get_ydim(...)";
		throw runtime_error(err_str);
	}

	return rg_ynum;
}

double gctl::regular_grid::get_xmin() const
{
	if (!initialized)
	{
		std::string err_str = "The object is not initialized. From void regular_grid::get_xmin(...)";
		throw runtime_error(err_str);
	}

	return rg_xmin;
}

double gctl::regular_grid::get_ymin() const
{
	if (!initialized)
	{
		std::string err_str = "The object is not initialized. From void regular_grid::get_ymin(...)";
		throw runtime_error(err_str);
	}

	return rg_ymin;
}

double gctl::regular_grid::get_dx() const
{
	if (!initialized)
	{
		std::string err_str = "The object is not initialized. From void regular_grid::get_dx(...)";
		throw runtime_error(err_str);
	}

	return rg_dx;
}

double gctl::regular_grid::get_dy() const
{
	if (!initialized)
	{
		std::string err_str = "The object is not initialized. From void regular_grid::get_dy(...)";
		throw runtime_error(err_str);
	}

	return rg_dy;
}

#ifdef GCTL_NETCDF

void gctl::regular_grid::load_netcdf_grid(std::string filename, mesh_data_type_e d_type, std::string xname, 
	std::string yname, std::string dataname)
{
	gctl::array<double> in_arr;
	gctl::array<double> in_x;
	gctl::array<double> in_y;

	// 如果微初始化 则初始化
	if (!initialized)
	{
		// read a data array
		read_netcdf_axis(filename, in_x, xname);
		read_netcdf_axis(filename, in_y, yname);
		read_netcdf_grid(filename, in_arr, xname, yname, dataname);

		int xnum = in_x.size();
		int ynum = in_y.size();
		double xmin = GCTL_BDL_MAX, ymin = GCTL_BDL_MAX;
		double xmax = GCTL_BDL_MIN, ymax = GCTL_BDL_MIN;

		for (int i = 0; i < xnum; i++)
		{
			xmin = GCTL_MIN(xmin, in_x[i]);
			xmax = GCTL_MAX(xmax, in_x[i]);
		}

		for (int i = 0; i < ynum; i++)
		{
			ymin = GCTL_MIN(ymin, in_y[i]);
			ymax = GCTL_MAX(ymax, in_y[i]);
		}

		double dx = (xmax - xmin)/(xnum - 1);
		double dy = (ymax - ymin)/(ynum - 1);

		meshdata *data_ptr;
		if (d_type == NodeData)
		{
			init(filename, "Imported from a .nc file.", xnum, ynum, xmin, ymin, dx, dy);
			data_ptr = meshdata_double::create(dataname, NodeData, xnum*ynum, true, 0.0);
		}
		else if (d_type == ElemData)
		{
			init(filename, "Imported from a .nc file.", xnum+1, ynum+1, xmin-0.5*dx, ymin-0.5*dy, dx, dy);
			data_ptr = meshdata_double::create(dataname, ElemData, xnum*ynum, true, 0.0);
		}

		saved_data.push_back(data_ptr);

		// get data pointer
		gctl::array<double> *val_ptr = (gctl::array<double>*) data_ptr->get_datval();
		for (int i = 0; i < val_ptr->size(); i++)
		{
			val_ptr->at(i) = in_arr[i];
		}

		initialized = true;
	}
	else
	{
		// read a data array
		read_netcdf_grid(filename, in_arr, xname, yname, dataname);

		std::string err_str;
		meshdata *data_ptr;
		if (d_type == NodeData)
		{
			if (in_arr.size() != rg_xnum*rg_ynum)
			{
				err_str = "Input data size is not compatible. From void regular_grid::load_netcdf(...)";
				throw runtime_error(err_str);
			}

			data_ptr = meshdata_double::create(dataname, d_type, rg_xnum*rg_ynum, true, 0.0);
		}
		else if (d_type == ElemData)
		{
			if (in_arr.size() != (rg_xnum-1)*(rg_ynum-1))
			{
				err_str = "Input data size is not compatible. From void regular_grid::load_netcdf(...)";
				throw runtime_error(err_str);
			}

			data_ptr = meshdata_double::create(dataname, d_type, (rg_xnum-1)*(rg_ynum-1), true, 0.0);
		}

		saved_data.push_back(data_ptr);

		// get data pointer
		gctl::array<double> *val_ptr = (gctl::array<double>*) data_ptr->get_datval();
		for (int i = 0; i < val_ptr->size(); i++)
		{
			val_ptr->at(i) = in_arr[i];
		}
	}
	return;
}

void gctl::regular_grid::save_netcdf_grid(std::string filename, mesh_data_type_e d_type)
{
	if (!initialized)
	{
		std::string err_str = "The object is not initialized. From void regular_grid::save_netcdf_grid(...)";
		throw runtime_error(err_str);
	}

	meshdata *curr_data = nullptr;

	bool if_append = false;
	for (iter = saved_data.begin(); iter != saved_data.end(); ++iter)
	{
		curr_data = *iter;
		if (curr_data->get_dattype() == d_type && d_type == NodeData && 
			curr_data->get_valtype() == Double && curr_data->get_output())
		{
			array<double>* data_ptr = (array<double>*)curr_data->get_datval();

			if (!if_append)
			{
				gctl::save_netcdf_grid(filename, *data_ptr, rg_xnum, rg_ynum, rg_xmin, 
					rg_dx, rg_ymin, rg_dy, "x", "y", curr_data->get_datname());
				if_append = true;
			}
			else
			{
				gctl::append_netcdf_grid(filename, *data_ptr, "x", "y", curr_data->get_datname());
			}
		}
		else if (curr_data->get_dattype() == d_type && d_type == ElemData && 
			curr_data->get_valtype() == Double && curr_data->get_output())
		{
			array<double>* data_ptr = (array<double>*)curr_data->get_datval();

			if (!if_append)
			{
				gctl::save_netcdf_grid(filename, *data_ptr, rg_xnum-1, rg_ynum-1, rg_xmin+0.5*rg_dx, 
					rg_dx, rg_ymin+0.5*rg_dy, rg_dy, "x", "y", curr_data->get_datname());
				if_append = true;
			}
			else
			{
				gctl::append_netcdf_grid(filename, *data_ptr, "x", "y", curr_data->get_datname());
			}
		}
	}

	return;
}

void gctl::regular_grid::save_netcdf_grid(std::string filename, std::string datname)
{
	meshdata *curr_data = get_data(datname);

	std::string err_str;
	if (!curr_data->get_output())
	{
		err_str = "Outputting of the data is disabled. From void regular_grid::save_netcdf_grid(...)";
		throw runtime_error(err_str);
	}

	if (curr_data->get_valtype() != Double)
	{
		err_str = "The function can not be applied to current data type. From void regular_grid::save_netcdf_grid(...)";
		throw runtime_error(err_str);
	}

	if (curr_data->get_valtype() == Double && curr_data->get_dattype() == NodeData)
	{
		array<double>* data_ptr = (array<double>*)curr_data->get_datval();
		gctl::save_netcdf_grid(filename, *data_ptr, rg_xnum, rg_ynum, rg_xmin, 
			rg_dx, rg_ymin, rg_dy, "x", "y", curr_data->get_datname());
	}
	else if (curr_data->get_valtype() == Double && curr_data->get_dattype() == ElemData)
	{
		array<double>* data_ptr = (array<double>*)curr_data->get_datval();
		gctl::save_netcdf_grid(filename, *data_ptr, rg_xnum-1, rg_ynum-1, rg_xmin+0.5*rg_dx, 
			rg_dx, rg_ymin+0.5*rg_dy, rg_dy, "x", "y", curr_data->get_datname());
	}

	return;
}

#endif // GCTL_NETCDF

void gctl::regular_grid::load_surfer_grid(std::string filename, std::string datname, mesh_data_type_e d_type, surfer_file_type_e grid_type)
{
	gctl::array<double> in_arr;

	int xnum, ynum;
	double xmin, xmax, ymin, ymax, zmin, zmax, dx, dy, blank_val;

	// 如果微初始化 则初始化
	if (!initialized)
	{
		if (grid_type == Surfer6Text)
		{
			gctl::read_surfer6_grid(filename, in_arr, xnum, ynum, xmin, xmax, ymin, ymax, zmin, 
				zmax, gctl::Surfer6Text);

			dx = (xmax-xmin)/(xnum-1);
			dy = (ymax-ymin)/(ynum-1);
		}
		else if (grid_type == Surfer6Binary)
		{
			gctl::read_surfer6_grid(filename, in_arr, xnum, ynum, xmin, xmax, ymin, ymax, zmin, 
				zmax, gctl::Surfer6Binary);

			dx = (xmax-xmin)/(xnum-1);
			dy = (ymax-ymin)/(ynum-1);
		}
		else if (grid_type == Surfer7Grid)
		{
			gctl::read_surfer7_grid(filename, in_arr, xnum, ynum, xmin, ymin, dx, dy, zmin, zmax, blank_val);
		}

		meshdata *data_ptr;
		if (d_type == NodeData)
		{
			init(filename, "Imported from a surfer grid.", xnum, ynum, xmin, ymin, dx, dy);
			data_ptr = meshdata_double::create(datname, NodeData, xnum*ynum, true, 0.0);
		}
		else if (d_type == ElemData)
		{
			init(filename, "Imported from a surfer grid.", xnum+1, ynum+1, xmin-0.5*dx, ymin-0.5*dy, dx, dy);
			data_ptr = meshdata_double::create(datname, ElemData, xnum*ynum, true, 0.0);
		}

		saved_data.push_back(data_ptr);

		// get data pointer
		gctl::array<double> *val_ptr = (gctl::array<double>*) data_ptr->get_datval();
		for (int i = 0; i < val_ptr->size(); i++)
		{
			val_ptr->at(i) = in_arr[i];
		}

		initialized = true;
	}
	else
	{
		if (grid_type == Surfer6Text)
		{
			gctl::read_surfer6_grid(filename, in_arr, xnum, ynum, xmin, xmax, ymin, ymax, zmin, 
				zmax, gctl::Surfer6Text);
		}
		else if (grid_type == Surfer6Binary)
		{
			gctl::read_surfer6_grid(filename, in_arr, xnum, ynum, xmin, xmax, ymin, ymax, zmin, 
				zmax, gctl::Surfer6Binary);
		}
		else if (grid_type == Surfer7Grid)
		{
			gctl::read_surfer7_grid(filename, in_arr, xnum, ynum, xmin, ymin, dx, dy, zmin, zmax, blank_val);
		}

		std::string err_str;
		meshdata *data_ptr;
		if (d_type == NodeData)
		{
			if (in_arr.size() != rg_xnum*rg_ynum)
			{
				err_str = "Input data size is not compatible. From void regular_grid::load_netcdf(...)";
				throw runtime_error(err_str);
			}

			data_ptr = meshdata_double::create(datname, d_type, rg_xnum*rg_ynum, true, 0.0);
		}
		else if (d_type == ElemData)
		{
			if (in_arr.size() != (rg_xnum-1)*(rg_ynum-1))
			{
				err_str = "Input data size is not compatible. From void regular_grid::load_netcdf(...)";
				throw runtime_error(err_str);
			}

			data_ptr = meshdata_double::create(datname, d_type, (rg_xnum-1)*(rg_ynum-1), true, 0.0);
		}

		saved_data.push_back(data_ptr);

		// get data pointer
		gctl::array<double> *val_ptr = (gctl::array<double>*) data_ptr->get_datval();
		for (int i = 0; i < val_ptr->size(); i++)
		{
			val_ptr->at(i) = in_arr[i];
		}
	}
	return;
}

void gctl::regular_grid::save_surfer_grid(std::string filename, std::string datname, surfer_file_type_e grid_type)
{
	meshdata *curr_data = curr_data = get_data(datname);

	std::string err_str;
	if (!curr_data->get_output())
	{
		err_str = "Outputting of the data is disabled. From void regular_grid::save_surfer_grid(...)";
		throw runtime_error(err_str);
	}

	if (curr_data->get_valtype() != Double)
	{
		err_str = "The function can not be applied to current data type. From void regular_grid::save_surfer_grid(...)";
		throw runtime_error(err_str);
	}

	gctl::array<double> *data_ptr = (gctl::array<double>*) curr_data->get_datval();
	if (curr_data->get_dattype() == NodeData)
	{
		if (grid_type == Surfer6Text)
		{
			gctl::save_surfer6_grid(filename, *data_ptr, rg_xnum, rg_ynum, 
				rg_xmin, rg_xmin+(rg_xnum-1)*rg_dx, rg_ymin, rg_ymin+(rg_ynum-1)*rg_dy, 
				NAN, NAN, gctl::Surfer6Text);
		}
		else if (grid_type == Surfer6Binary)
		{
			gctl::save_surfer6_grid(filename, *data_ptr, rg_xnum, rg_ynum, 
				rg_xmin, rg_xmin+(rg_xnum-1)*rg_dx, rg_ymin, rg_ymin+(rg_ynum-1)*rg_dy, 
				NAN, NAN, gctl::Surfer6Binary);
		}
		else if (grid_type == Surfer7Grid)
		{
			gctl::save_surfer6_grid(filename, *data_ptr, rg_xnum, rg_ynum, 
				rg_xmin, rg_ymin, rg_dx, rg_dy);
		}
	}
	else if (curr_data->get_dattype() == ElemData)
	{
		if (grid_type == Surfer6Text)
		{
			gctl::save_surfer6_grid(filename, *data_ptr, rg_xnum-1, rg_ynum-1, 
				rg_xmin+0.5*rg_dx, rg_xmin+0.5*rg_dx+(rg_xnum-2)*rg_dx, 
				rg_ymin+0.5*rg_dy, rg_ymin+0.5*rg_dy+(rg_ynum-2)*rg_dy, 
				NAN, NAN, gctl::Surfer6Text);
		}
		else if (grid_type == Surfer6Binary)
		{
			gctl::save_surfer6_grid(filename, *data_ptr, rg_xnum-1, rg_ynum-1, 
				rg_xmin+0.5*rg_dx, rg_xmin+0.5*rg_dx+(rg_xnum-2)*rg_dx, 
				rg_ymin+0.5*rg_dy, rg_ymin+0.5*rg_dy+(rg_ynum-2)*rg_dy, 
				NAN, NAN, gctl::Surfer6Binary);
		}
		else if (grid_type == Surfer7Grid)
		{
			gctl::save_surfer6_grid(filename, *data_ptr, rg_xnum-1, rg_ynum-1, 
				rg_xmin+0.5*rg_dx, rg_ymin+0.5*rg_dy, rg_dx, rg_dy);
		}
	}

	return;
}

void gctl::regular_grid::save_gmsh(std::string filename, index_packed_e packed)
{
	std::ofstream outfile;
	gctl::open_outfile(outfile, filename, ".msh");
	gctl::save2gmsh(outfile, elements, nodes, packed);
	outfile.close();
	return;
}

void gctl::regular_grid::save_gmsh(std::string filename, mesh_data_type_e d_type, output_type_e out_mode, index_packed_e packed)
{
	base_mesh::save_gmsh(filename, d_type, out_mode, packed);
	return;
}

void gctl::regular_grid::save_gmsh(std::string filename, std::string datname, output_type_e out_mode, index_packed_e packed)
{
	base_mesh::save_gmsh(filename, datname, out_mode, packed);
	return;
}

void gctl::regular_grid::load_data_cloud(const array<point2dc> &in_posi, const array<double> &in_val, 
	double search_xlen, double search_ylen, double search_deg, std::string datname, mesh_data_type_e d_type)
{
	gctl::array<gctl::point2dc> grid_points;

	std::string err_str;
	if (in_posi.size() != in_val.size())
	{
		err_str = "Size of the input arrays are different. From gctl::regular_grid::load_data_cloud(...)";
		throw runtime_error(err_str);
	}
	// create a new data object
	meshdata *data_ptr = add_data(datname, d_type, true, Double);

	array<double> *inte_data = nullptr;
	if (d_type == NodeData)
	{
		grid_points.resize(rg_xnum*rg_ynum);
		for (int j = 0; j < rg_ynum; j++)
		{
			for (int i = 0; i < rg_xnum; i++)
			{
				grid_points[i+j*rg_xnum].x = rg_xmin+i*rg_dx;
				grid_points[i+j*rg_xnum].y = rg_ymin+j*rg_dy;
			}
		}
	}
	else if (d_type == ElemData)
	{
		grid_points.resize((rg_xnum-1)*(rg_ynum-1));
		for (int j = 0; j < rg_ynum-1; j++)
		{
			for (int i = 0; i < rg_xnum-1; i++)
			{
				grid_points[i+j*(rg_xnum-1)].x = rg_xmin+i*rg_dx+0.5*rg_dx;
				grid_points[i+j*(rg_xnum-1)].y = rg_ymin+j*rg_dy+0.5*rg_dy;
			}
		}
	}

	inte_data = geometry2d::p2p_dist_2d(grid_points.get(), grid_points.size(), in_posi.get(), 
		in_val.get(), in_posi.size(), search_xlen, search_ylen, search_deg);

	array<double> *val_ptr = (array<double>*) data_ptr->get_datval();
	for (int i = 0; i < val_ptr->size(); i++)
	{
		val_ptr->at(i) = inte_data->at(i);
	}

	delete inte_data;
	return;
}

void gctl::regular_grid::extract_points(std::string datname, const array<point2dc> &in_posi, 
	array<double> &out_val)
{
	std::string err_str;
	if (in_posi.size() == 0)
	{
		err_str = "The input array is empty. From gctl::regular_grid::extract_points(...)";
		throw runtime_error(err_str);
	}

	meshdata *data_ptr = get_data(datname);
	mesh_data_value_e value_type = data_ptr->get_valtype();
	mesh_data_type_e data_type = data_ptr->get_dattype();
	if(value_type != Double)
	{
		err_str = "Wrong data's value type. From gctl::regular_grid::extract_points(...)";
		throw runtime_error(err_str);
	}

	array<double> *val_ptr = (array<double>*) data_ptr->get_datval();
	out_val.resize(in_posi.size());

	int tmp_i, tmp_j, in_xnum, in_ynum;
	if (data_type == NodeData)
	{
		in_xnum = rg_xnum;
		in_ynum = rg_ynum;
		for (int i = 0; i < in_posi.size(); i++)
		{
			tmp_i = floor((in_posi[i].x - rg_xmin)/rg_dx);
			tmp_j = floor((in_posi[i].y - rg_ymin)/rg_dy);

			if (tmp_i > in_xnum-2 || tmp_i < 0 || tmp_j > in_ynum-2 || tmp_j < 0)
			{
				out_val.at(i) = NAN;
			}

			if (!std::isnan(val_ptr->at(tmp_i + tmp_j*in_xnum)) && 
				!std::isnan(val_ptr->at(tmp_i + 1 + tmp_j*in_xnum)) && 
				!std::isnan(val_ptr->at(tmp_i + 1 + (tmp_j+1)*in_xnum)) && 
				!std::isnan(val_ptr->at(tmp_i + (tmp_j+1)*in_xnum)))
			{
				out_val.at(i) = rect_interpolate(rg_xmin + tmp_i*rg_dx, rg_ymin + tmp_j*rg_dy, 
					rg_dx, rg_dy, in_posi[i].x, in_posi[i].y, val_ptr->at(tmp_i + tmp_j*in_xnum), 
					val_ptr->at(tmp_i + 1 + tmp_j*in_xnum), val_ptr->at(tmp_i + 1 + (tmp_j+1)*in_xnum), 
					val_ptr->at(tmp_i + (tmp_j+1)*in_xnum));
			}
			else out_val.at(i) = NAN;
		}
	}
	else
	{
		in_xnum = rg_xnum-1;
		in_ynum = rg_ynum-1;
		for (int i = 0; i < in_posi.size(); i++)
		{
			tmp_i = floor((in_posi[i].x - rg_xmin - 0.5*rg_dx)/rg_dx);
			tmp_j = floor((in_posi[i].y - rg_ymin - 0.5*rg_dy)/rg_dy);

			if (tmp_i > in_xnum-2 || tmp_i < 0 || tmp_j > in_ynum-2 || tmp_j < 0)
			{
				out_val.at(i) = NAN;
			}

			if (!std::isnan(val_ptr->at(tmp_i + tmp_j*in_xnum)) && 
				!std::isnan(val_ptr->at(tmp_i + 1 + tmp_j*in_xnum)) && 
				!std::isnan(val_ptr->at(tmp_i + 1 + (tmp_j+1)*in_xnum)) && 
				!std::isnan(val_ptr->at(tmp_i + (tmp_j+1)*in_xnum)))
			{
				out_val.at(i) = rect_interpolate(rg_xmin + 0.5*rg_dx + tmp_i*rg_dx, rg_ymin + 0.5*rg_dy + tmp_j*rg_dy, 
					rg_dx, rg_dy, in_posi[i].x, in_posi[i].y, val_ptr->at(tmp_i + tmp_j*in_xnum), 
					val_ptr->at(tmp_i + 1 + tmp_j*in_xnum), val_ptr->at(tmp_i + 1 + (tmp_j+1)*in_xnum), 
					val_ptr->at(tmp_i + (tmp_j+1)*in_xnum));
			}
			else out_val.at(i) = NAN;
		}
	}
	return;
}

void gctl::regular_grid::extract_profile(std::string datname, const point2dc &start_p, 
	const point2dc &end_p, int size_p, array<point2dc> &out_posi, array<double> &out_val)
{
	std::string err_str;
	if (start_p == end_p)
	{
		err_str = "The start and end points can's be the same. From gctl::regular_grid::extract_profile(...)";
		throw runtime_error(err_str);
	}

	if (size_p <= 0)
	{
		err_str = "Invalid profile size. From gctl::regular_grid::extract_profile(...)";
		throw runtime_error(err_str);
	}

	linespace(start_p, end_p, size_p, out_posi);
	extract_points(datname, out_posi, out_val);
	return;
}

void gctl::regular_grid::gradient(std::string datname, std::string gradname, gradient_type_e d_type, int order)
{
	meshdata *data_ptr = get_data(datname);
	mesh_data_value_e value_type = data_ptr->get_valtype();
	mesh_data_type_e data_type = data_ptr->get_dattype();

	std::string err_str;
	if(value_type != Double)
	{
		err_str = "Wrong data's value type. From gctl::regular_grid::gradient(...)";
		throw runtime_error(err_str);
	}

	if (d_type != Dx && d_type != Dy)
	{
		err_str = "The calculation type must be Dx or Dy. From gctl::regular_grid::gradient(...)";
		throw runtime_error(err_str);
	}

	// 检查是否存在与梯度数据同名的数据 若有则检查是否需要重新建立数据
	meshdata *new_data_ptr = nullptr;
	mesh_data_value_e value_type2;
	mesh_data_type_e data_type2;
	if (saved(gradname))
	{
		new_data_ptr = get_data(gradname);
		value_type2 = new_data_ptr->get_valtype();
		data_type2  = new_data_ptr->get_dattype();
		if (value_type2 != value_type || data_type2 != data_type)
		{
			// 删除原有数据
			remove_data(gradname);
			// 新建数据
			new_data_ptr = add_data(gradname, data_type, true, value_type);
		}
	}
	else new_data_ptr = add_data(gradname, data_type, true, value_type);

	array<double> *inval_ptr = (array<double>*) data_ptr->get_datval();
	array<double> *outval_ptr = (array<double>*) new_data_ptr->get_datval();
	if (data_type == NodeData && d_type == Dx)
	{
		difference_2d(*inval_ptr, *outval_ptr, rg_ynum, rg_xnum, rg_dx, Dx, order);
	}
	else if (data_type == NodeData && d_type == Dy)
	{
		difference_2d(*inval_ptr, *outval_ptr, rg_ynum, rg_xnum, rg_dy, Dy, order);
	}
	else if (data_type == ElemData && d_type == Dx)
	{
		difference_2d(*inval_ptr, *outval_ptr, rg_ynum-1, rg_xnum-1, rg_dx, Dx, order);
	}
	else
	{
		difference_2d(*inval_ptr, *outval_ptr, rg_ynum-1, rg_xnum-1, rg_dy, Dy, order);
	}
	return;
}

#ifdef GCTL_WAVELIB

void gctl::regular_grid::wavelet(std::string datname, std::string wavename, int order, bool summary)
{
	meshdata *data_ptr = get_data(datname);
	mesh_data_value_e value_type = data_ptr->get_valtype();
	mesh_data_type_e data_type = data_ptr->get_dattype();

	std::string err_str;
	if(value_type != Double)
	{
		err_str = "Wrong data's value type. From gctl::regular_grid::wavelet(...)";
		throw runtime_error(err_str);
	}

	if (order <= 0)
	{
		err_str = "The decomposition level must be greater than zero. From gctl::regular_grid::wavelet(...)";
		throw runtime_error(err_str);
	}

	array<double> *inval_ptr = (array<double>*) data_ptr->get_datval();

	wave_object obj = nullptr;
	wt2_object wt = nullptr;
	obj = wave_init(wavename.c_str());
	if (data_type == NodeData) wt = wt2_init(obj, "dwt", rg_ynum, rg_xnum, order);
	else wt = wt2_init(obj, "dwt", rg_ynum-1, rg_xnum-1, order);

	double *wavecoeffs = nullptr;
	wavecoeffs = dwt2(wt, inval_ptr->get());
	gctl::array<double> coeff_copy(wt->outlength);

	for (int i = 0; i < wt->outlength; i++)
	{
		coeff_copy[i] = wavecoeffs[i];
		wavecoeffs[i] = 0.0;
	}

	int ir = wt->dimensions[0], ic = wt->dimensions[1];
	for (int i = 0; i < ir*ic; i++)
	{
		wavecoeffs[i] = coeff_copy[i];
	}

	data_ptr = add_data(datname+"_A"+std::to_string(wt->J), data_type, true, value_type);
	array<double> *outval_ptr = (array<double>*) data_ptr->get_datval();
	idwt2(wt, wavecoeffs, outval_ptr->get());

	int coeff_idx;
	for (int j = 0; j < wt->J; j++)
	{
		for (int i = 0; i < wt->outlength; i++)
		{
			wavecoeffs[i] = 0.0;
		}

		ir = wt->dimensions[2*j];
		ic = wt->dimensions[2*j+1];
		coeff_idx = wt->coeffaccess[3*j+1];
		for (int i = 0; i < 3*ir*ic; i++)
		{
			wavecoeffs[i+coeff_idx] = coeff_copy[i+coeff_idx];
		}

		data_ptr = add_data(datname+"_D"+std::to_string(wt->J - j), data_type, true, value_type);
		outval_ptr = (array<double>*) data_ptr->get_datval();
		idwt2(wt, wavecoeffs, outval_ptr->get());
	}

	if (summary) wt2_summary(wt);

	if (obj != nullptr) wave_free(obj);
	if (wt != nullptr) wt2_free(wt);
	if (wavecoeffs != nullptr) free(wavecoeffs);
	return;
}

#endif // GCTL_WAVELIB
